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Abstract. The steady-state distributions and dynamical behaviour of Zero Range 
Processes with hopping rates which are non-monotonic functions of the site occupation 
are studied. We consider two classes of non-monotonic hopping rates. The first results 
in a condensed phase containing a large (but subextensive) number of mesocondensates 
each containing a subextensive number of particles. The second results in a condensed 
phase containing a finite number of extensive condensates. We study the scaling 
behaviour of the peak in the distribution function corresponding to the condensates 
in both cases. In studying the dynamics of the condensate we identify two timescales: 
one for creation, the other for evaporation of condensates at a given site. The scaling 
behaviour of these timescales is studied within the Arrhenius law approach and by 
numerical simulations. 
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1. Introduction 

Real-space condensation is a phenomenon which occurs in a variety of physical systems 
such as jamming in traffic flow, granular clustering, wealth condensation and hub 
formation in complex networks [1, 2]. In a condensation process a finite fraction of 
the microscopic constituents aggregate in space. A simple (minimal) model which has 
commonly been used to describe this phenomenon in recent years, is the Zero-Range 
Process (ZRP)[3]. The model comprises particles distributed over the sites of a lattice, 
with stochastic dynamics allowing particles to hop between sites. The hopping rate is 
just a function of the occupation, n, of the site from which a hop takes place, hence the 
name ZRP. For recent reviews of the properties of the ZRP see [4, 5]. 

An appealing feature of the ZRP is that its steady-state distribution is exactly 
calculable in terms of the hopping rates. It has been shown that for a class of hopping 
rates which are monotonic and decrease suitably slowly with n, the model exhibits 
a condensation transition at some critical value of the particle density (the average 
number of particles per lattice site). Below the critical density the system is in the fluid 
phase where all sites are equivalent and have a characteristic occupation given by the 
average density. Above the critical density the excess particles — a finite fraction of the 
total number of particles — condense onto a single site. The supercritical system is thus 
composed of a critical fluid coexisting with a condensate. 

A choice of hopping rate commonly used to study condensation is 

u(n) = 1 + - for n > , (1) 
u(0) = . 

For b < 2 it has been shown that there is no condensation at any density, while for b > 2 
condensation takes place above a critical density [4, 5, 6, 7] 

* = f^- (2) 

Our interest in this work is to consider dynamics which could lead to the formation 
of more than one condensate. Previous studies have shown that multiple condensates can 
exist in models with non-conserving dynamics where in addition to the above hopping 
rate (1), local creation and annihilation processes are introduced [8, 9]. Here we consider 
ZRP models with conserving dynamics (i.e. where particles are neither created nor 
destroyed) which generate multiple condensates. A simple mechanism for this is to 
choose hopping rates which are no longer monotonic but increase for large n. This 
effectively introduces a soft cut-off in the occupation number of a site. Then, the excess 
particles are forced to be shared over multiple condensate sites. We study both steady- 
state properties and the dynamical features corresponding to creation and evaporation 
of condensates in such a model. 

The paper is organised as follows. The model and different dynamics we consider 
are defined and general steady state properties are reviewed in Section 2. In Section 3 
we analyse the properties of the steady state in the case of an algebraically increasing 
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cut-off in the hop rates. In Section 4 we analyse the properties of the steady state in 
the case of an exponentially increasing cut-off in the hop rates. In particular the scaling 
of the condensate properties are compared with numerical simulations. The dynamics 
of the condensates are studied in Section 5. Finally, some conclusions are drawn in 
Section 6. 

2. Model definition 

We consider a lattice of L sites with a total of N particles and average density p = N/ L. 
Each occupation of a lattice site / is given by n\. One choice of hopping rate we consider 
is given by 

b / n\ k 

u(n) = 1 + - + c - for n > . (3) 
n \L J 

With this choice, illustrated in Figure 1, the hopping rate from a site containing 
n = O(L) particles is enhanced over the case (I). The hopping rate is non-monotonic: 
initially it decreases as 1/n, it has a minimum at n ~ j^k/k+i an j ^ ^hen increases 
algebraically. We refer to the choice (3) as algebraic hopping rates. We note that the 
hopping rates (3) depend on the system size; another ZRP with size-dependent hopping 
rates that induce condensation has been studied in [14]. 

For b < 2, where condensation does not occur, the additional term in (3) compared 
to (1) does not affect the distribution of occupations. On the other hand for b > 2, we 
expect condensation to take place at the same critical density as for hopping rates (3). 
The reason is that the last term in (3) is negligibly small for nonextensive occupations, 
therefore does not affect the fluid phase. Hence, the condensation transition will be at 
the same critical density p c . However the last term in (3) does affect the condensed 
phase since the last term becomes significant for macroscopic n and thus would change 
the occupation of any condensate. We shall show that, in fact, in the condensed phase 
the system exhibits a large number of mesocondensates whose occupation increases 
sublinearly with L. The number of mesocondensates is also sublinear in L but the total 
occupation of all condensates is extensive, accounting for the total number of excess 
particles. 

Another choice of hopping rates that we consider is 

u{n) = 1 + — + exp(n — aL) for n > , (4) 

where a > 0. The hopping rate is again non-monotonic as illustrated in Figure 1: 
initially it decreases as 1/n, it has a minimum at n ~ aL and it then increases 
exponentially. We refer to the choice (4) as exponential hopping rates. Here the increase 
in the hopping rate with n at n = aL is sharper than (3). In this case the effective cut-off 
is at an extensive occupation and, when present, condensates contain up to aL particles. 
Thus, at sufficiently high density a finite number of condensates each with occupancy 
< aL may be present. 
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Figure 1. The algebraic hopping rate u(n) given by (3) for b = 4, L = 1000 and 
different values of k and c. The full line represents c — where u(n) is a monotonically 
decreasing function of n whereas the dashed and dotted lines represent c = 1 and 
k = 1, 10 respectively where u(n) has a minimum. The exponential hopping rate given 
by (4) is represented by the dashed-dotted line . 

In a numerical simulation of the dynamics, the hopping rates are implemented in 
a random sequential updating scheme by first selecting a site at random. Then one 
should generate a random number uniformly distributed between and the maximum 
possible hopping rate and carry out a move if the random number is less than or equal 
to the rate. If the maximal rate is u(N), this rate may not be realised in practice since 
it would correspond to all particles residing in a single site. Therefore in order to speed 
up the simulation we choose a cut-off in the hopping rate so that the dynamics is not 
affected. 

When a hopping process is carried out the destination site must be chosen. In 
this work we restrict ourselves to one-dimensional totally asymmetric hopping where 
the sites are arranged on a ring and particles hop from site % to the next site % + 1. In 
general, the steady-state distribution for the ZRP does not depend on the connectivity 
although the dynamical properties may do. 

2.1. Steady-State Properties 

The steady state of a ZRP is known exactly and has a factorised form: 

P( ni , ■■■,n L ) = Z-\L, N) n f(rn)6(jr nj - N) (5) 

i=i j=i 
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Figure 2. A logarithmic plot of the single site distribution p(n) for L = 1600, 
N = 6400 calculated from the exact expression (7). The hopping rate u{n) is given by 
the algebraic case (3) with b = 4, c = 1 and k = 10. 



where the normalisation Z(L,N), equivalent to a canonical partition function for N 
particles and L sites, is defined as 



Z(L,N) 



L 

n 



N 

E /fa) 

71, =0 



TV) 



(6) 



where the delta symbol ensures that only configurations with precisely N particles 
contribute. The single site distribution, which is the probability that in the steady 
state a given site contains n particles, is given by 

Z(L-l,N-n) 



Pin) = f(n)'- 



(7) 



Z(L,N) 

To compute Z(L, N) and p(n) exactly for small systems the following recursion is useful 



A? 



Z(L,N) = J2f(n)Z(L-l,N-n). 



(8) 



n=0 



We will use this formulation (namely, the Canonical Ensemble formulation in which 
the total number of particles is fixed) to provide numerical results. This data will be 
compared in Sections 3.2 and 4.2 with the analytical predictions of Sections 3.1 and 4.1 
that are made within the grand canonical ensemble which we now describe. 

To describe analytically the properties of p(n) it is convenient to work within the 
grand canonical ensemble (where the total number of particles is allowed to fluctuate 
and one introduces a fugacity z to control the number of particles). The steady-state 
distribution becomes 

L 

P(n ir --,n L ) = Hp(ni) , (9) 



where 



p(n) = Az n f(n) , 



(10) 
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and 

n i 

/(») = n -n ( n ) 

mil < m ) 

The normalization constant A in (10) is chosen to ensure that sum of p(n) is equal to 1. 
The fugacity z is determined by the condition that the mean total number of particles 
in the system is N. 

In the commonly studied case with u{n) given by (1), the probability distribution 
at large n becomes 

p(n) ~ A Z — . (12) 

The condensation transition is easy to understand. For densities below the critical 
density, p c , the fugacity satisfies z < 1 and the probability decays exponentially for large 
n. This distribution describes what is termed the fluid phase. As the critical density is 
approached, z approaches 1 and the decay of p{n) becomes algebraic, corresponding to 
a critical fluid. Above the critical density an extra piece of p(n) emerges representing 
the condensate [10, 11]. Therefore the condensate co-exists with the critical fluid which 
has density p c . The number of particles in the condensate is given by L(p — p c ). 



3. Analysis of condensation for algebraically increasing u(n) 

In this Section we study the characteristic features of the single site occupation 
distribution p(n) for the case of algebraically increasing hopping rate (3). In the 
condensed phase we expect this distribution to have the general form given in Figure 2. 
The distribution decays algebraically for intermediate n and there is a peak in the 
distribution for large n. The peak is associated with any condensates in the system. 
Our aim is to analyse the shape of this peak — its height, width and position. This 
information will determine the typical number of condensates present in the system and 
their typical occupancy. 



3.1. Grand Canonical Analysis 

To investigate the probability distribution in the condensed phase within the grand 
canonical ensemble we consider the effective potential $(n) defined through 

pin) = Ae-^ n) . (13) 

Using (10,11), $(n) is given by 

n 

<3>(n) = -n\n z + hiw(m) . (14) 

m=l 

Expanding lnu(m) for large m and approximating for n L the sum as an integral 
yields 

" c n k+1 

lnu(m) ~ - Hnn+ . (15) 

m=l K -\- I Li 
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Since we are discussing condensation we consider only hopping rates with b > 2 for 
which condensation is possible. In the condensed phase we expect the fugacity z — > 1 in 
the thermodynamic limit. Hence, one describes the fugacity as a function of the system 
size L, such that 

z = e h{L) , (16) 

where h(L) — > as L — > oo. 

Putting all this together we obtain, for large n and L, 



$(n) = —nh(L) + b Inn + 



n 



fc + 1 L fc ' 



which implies 



A 

- — exp 



n/i(L) 



c n 



fc+i 



;i7) 



;is) 



fc + l L fc 

where A is a normalization constant. 

We now identify the position of the peak which corresponds to a maximum of p(n) 
at n*. The condition 



$' W = -h(L) + - + = , 



gives 



Then /i(L) can be eliminated in (17) resulting in 

ck {n*) k+1 



+ b In n* — b 



(19) 



(20) 



(21) 



k + 1 L k 

Looking for a solution such that p(n*) scales as a power law of L implies that 
<3>(n*) = O(lnL). The leading order term in Equation (21) is the first term, and therefore 
n *^ L k/(k + i) [lnL ]V(fc+i) ; (22) 

which implies that the leading order of h(L) is 



h(L) 



'\nL_y k+1 



The width of the peak at n* is given by 

-1/2 

An* = $"(n*) 



where 



$"(n) = ck 



TV 



k-1 



(23) 



(24) 



(25) 



Using the above expression for n*, and noting that the first term of this expression is 
the leading one for large L, we find 

An* ~ fc -i/2 L fe/(fc + i) [lnL] -(k-i)/2(k + i) (26) 
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Since An*/n* ~ l/llnL] 1 ^ 2 vanishes in the large L limit we expect the canonical and 
grand canonical ensembles to be equivalent. 

In order to calculate p(n*) we note that the weight, w, of the peak, that is the 
probability contained in the peak, is given by 

w~p( n *)An*. (27) 

In the condensed phase, the weight w represents the fraction of sites which are 
condensate sites. The fraction of particles contained within the condensates is then 
n*w. Since in the condensed phase this fraction is finite (i.e. p — p c ) we arrive at the 
condition 

p{n*)n*An* = 0(1) (28) 

yielding 

pin") ~ fcV2 L -2fc/(fe+l) [lnL ](fc-3)/2( fe+ l) (2Q) 

The weight of the peak, w thus scales as 

w ^L" fc /( fc+1 )[lnL] (1/(fc+1) . (30) 

The above analysis of the scaling of the peak in p(n) implies the following 
condensation behaviour. We first note that the size of the condensates, (22), scales 
sub-linearly with the system size and therefore these are termed mesocondensates. The 
typical number of mesocondensates is Lw which scales as up to logarithmic 

corrections. Therefore their number diverges sublinearly with L. 

We conclude by investigating the dip in p(n) to the left of n* denoted n min in 
Figure 2. As we discuss in Section 5 the dip probability is significant in determining 
the dynamics of the mesocondensates. Balancing the first two terms in the extremum 
condition (19) yields 

n mm ^j^-L k ^[\nL]- k ^ . (31) 

Then, p{n m i n ) is readily evaluated as 

p(n mm )^L- bk /^[\nL] bk/{k+1) . (32) 

In summary, the leading behaviour in the system size L of the shape of peak is 
given by equations (22,26,29) and the dip by (31,32). 

The analysis of this subsection has been carried out within the grand canonical 
ensemble. In the following subsection we will compare it with exact numerical results 
for finite systems obtained in the canonical ensemble. This will provide evidence that 
at least the leading behaviour for the scaling properties of the condensate peak are 
correctly given by (22,26,29). 
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Figure 3. A logarithmic plot of p(n) for L = 1000 calculated from the exact expression 
(7). The hopping rate u(n) is given by the algebraic case (3) with b = 3, p = 2 and 
various values of c and k. The dashed line corresponds to c = 0, the dotted line to 
c = 1, k — 1 and the full line to c = 1, k = 10 



3.2. Numerical results 

For u(n) given by (3) and for p > p c , we plot in Figure 3 the exact p(n) for L = 1000. 
This was computed by iterating the recursion relation (8) and using expression (7). 
As expected the probability distribution exhibits a power law decay representing the 
fluid and a peak corresponding to the condensates. In the c = case, there is a single 
condensate containing L[p — p c ] particles which results in a peak location that scales as 
n* ~ L. For b > 3, the second moment is finite and the peak has a width of An* ~ L 1 / 2 . 
For 2 < b < 3, the regime of anomalous fluctuations, the second moment diverges and 
the width of the peak scales as An* ~ L 1 /^ 1 ). For both regimes the peak is sharp in 
the sense that An* /n* — > 0. It is clearly seen in Figure 3 that the peak for the case c = 
is sharper and at larger n than for c > as expected from the above and (22,26,29). 

In Figure 4 we present some snapshots of the occupations of the sites in the steady 
state as obtained from numerical simulations. It can be seen that while for c = (where 
(3) reduces to (1)) there is a single condensate, for c > the excess density above the 
critical value is distributed among a number of mesocondensates. It is the distribution 
of the occupations of these mesocondensates that the peak in p{n) describes. 

Since the steady state properties of the distribution given in Section 3 were 
calculated in the grand canonical ensemble,we wish to compare the distribution obtained 
in the canonical ensemble with the one obtained in the grand canonical ensemble. In 
Figure 5 the probability distribution for the parameters L = 1600, 6 = 4, p = 4, c = 1 
and k = 5 was calculated numerically both in the grand-canonical ensemble (10) and in 
the canonical ensemble (7). It can be seen in Figure 5 that the distributions obtained 
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site i site i site i 



Figure 4. Site occupations in snapshots of the system obtained from numerical 
simulation for L = 1000, p = 2, in the case of algebraic hop rates with 6 = 3 and 
several values of c and k: from left to right for c = 1, k = 1, c = 1, k = 10 and c = 0. 




n 



Figure 5. A comparison of the canonical probability (7) and the grand canonical 
probability (10), represented by the full line and (o) respectively. The distributions 
were calculated numerically for algebraic hopping (3) rates with L — 1600, b = 4, p = 4 
c = 1 and k = 5. 

form the two ensembles agree. 

In Figure 6 we check that the grand canonical analysis yields the leading scaling 
properties of the peak calculated for the canonical distribution in finite size systems. In 
Figure 6(a), p(n*) is plotted against L. A slope which agrees with the expected value of 
— 1.65 « —2k/(k + l) is seen, confirming the L _2fc /( fe+1 ) dependence (29). In Figure 6(b) 
n* is plotted against L. A slope which agrees with the expected value of 0.87 w kj (k + 1) 
is seen, confirming the L k ^ k+1 ^ dependence (22). The logarithmic corrections were not 
taken into account since the correction [lnL] 1 ^ 6 has negligible L dependence. 

In Figure 7 we check that the grand canonical analysis yields the leading scaling 
properties of the dip calculated for the canonical distribution in finite size systems. In 
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Figure 6. The scaling of the condensate peak as a function of the system size L in 
the case of algebraic hopping rates with b = 4, p = 4 and c = 1, k = 5. The probability 
p(n*) and location of the dip n* were calculated numerically for each system size L 
using (7) for the hopping rate (3). In (a) and (b) p{n*) and n* are plotted respectively 
as a a function of L and are compared to a linear fit of the data points. The data 
points for the various system sizes are compared to linear fit with a slope of —1.65 and 
0.87 in (a) and (b) respectively. 




L/ln(L) L/ln(L) 

Figure 7. The scaling of the dip as a function of the system size L in the case of 
algebraic hopping rates with b = 4, p = 4 and c = 1, k = 5. The probability p(n min ) 
and location of the dip n m i n were calculated numerically for each system size L using 
(7) for the hopping rate (3). In (a) and (b) p(n m i n ) and n m i„ are plotted respectively 
as a a function of L/\n(L) and are compared to a linear fit of the data points. The 
data points for the various system sizes are compare to linear slopes of —3.2 and 0.83 
in (a) and (b) respectively. 



Figure 7(a) the dip height p(n min ) is plotted against L/ln(L). A linear fit for the data 
with a slope of —3.2 « — bk/(k + l) is seen, confirming the [L/ \n(L)}~ bk ^ k+1 ^ dependence 
(32). In Figure 7(b) the dip location n min is plotted against Lj ln(L). A linear fit for the 
data with a slope of 0.83 w k/(k + l) is seen, confirming the [L/ ln(L)] fe ^ fc+1 - ) dependence 
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(31). 

We conclude that the steady state properties of the distribution given in Section 3 
computed for the grand canonical ensemble are able to describe the properties in the 
canonical ensemble. 



4. Analysis of condensation for exponentially increasing u(n) 

In this section we consider the distribution p(n) corresponding to hopping rates (4). As 
can be seen in Figure 1, the hopping rate u{n) is basically given by the commonly used 
form (1) but with a rapid increase at n — aL. In turn this generates a sharp cut-off in 
p{n) at n — aL. 

In the absence of the cut-off the condensate peak is at n = (p — p c )L. Thus, we 
expect this still to be the case as long as (p—p c ) < a. However, for densities above p c + a, 
and below p c + 2a, the density excess p — p c is distributed among two condensates each 
containing a maximum of aL particles. In general for densities p c +qa < p < p c +(q+l)a, 
where q is an integer, we would expect to have q + 1 condensates. In other words the 
number of condensates is 

9 + 1= T^l, (33) 

where \x] denotes the ceiling function which is the integer part of x plus one. Thus, 
in contrast to the algebraic case, this model exhibits a finite number of extensive 
condensates. 




Figure 8. The time-dependent occupation of the two condensates for exponential 
hopping rates with parameters L = 1000, a = 1, b = 4 and p = 2 calculated from 
numerical simulation. 
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For a density p = p c + qa with integer q, the condensed phase is composed of q 
condensates with aL particles each. The structure of the condensed phase for noninteger 
(p — Pc)/a is rather different. This is best illustrated by considering the case of two 
condensates: p c + a < p < p c + 2a. The total number of particles in the two condensate 
sites is (p — p c )L. Since the maximum occupation of a site is aL, the occupation of a 
condensate site fluctuates between aL and (p — p c — a)L. A numerical simulation was 
carried out for L = 1000, a = 1, b = 4 and p = 2 for which the expected number 
of condensates is \p — p c ] = 2. In Figure 8 the occupation of the two condensates is 
plotted as a function of time. Interestingly, the two condensates' occupations fluctuate 
in an anticorrelated fashion between aL and (p — p c — a)L around the mean occupation 
(p — p c )L/2. As discussed in [12], the condensates interact and the role of the larger 
condensate is shared among the two condensates. This results in a distribution with 
two peaks with one centred at n = aL corresponding to [(p — Pc)/a\ condensates and 
the second corresponds to a single condensate containing the remaining (p — p c — a)L 
particles. Here [^J is the floor function, defined as the highest integer less than or equal 
to x. The numerical results presented in section 4.2 confirm this behaviour. 

In the following subsection we analyze the scaling properties of the condensates 
within the grand canonical ensemble. The grand canonical analysis of the ZRP with 
a sharp cut-off in occupation, as is the case here, leads to a peak in the distribution 
p(n) at the cut-off for any density above p c [4] . On the other hand the considerations of 
the previous paragraph indicate that the behaviour of the model with fixed number of 
particles (Canonical ensemble) is more complex with a distribution depending on p. For 
example, as pointed out in the previous paragraph for p — p c < a the condensate peak 
is below the cut-off. We conclude that at densities for which (p — p c )/a is non-integer 
the grand canonical and canonical analyses do not yield the same scaling behaviour 
for the condensate. Thus we apply the grand canonical analysis for the case of integer 
(p — p c )/a where the condensate peak is at the cut-off point n = aL. In subsection 4.2 
we numerically check the equivalence between the two ensembles in this case. 



4-1. Grand Canonical Analysis 

Following the same strategy as in subsection 3.1, one finds that for hopping rate (4) 

$(n) = -nh(L) + 6 Inn + e n ~ aL , (34) 

which implies 



p{n) ~ ~ ex P nh(L) — e 

where A is a normalization constant. 

The extremum condition on is 

b 



n—aL 



$'( n ) = —h(L) + - + e n - aL = , 
n 



which implies 



$(n*) = -(n* - l)e n ~- aL + blnn* - b . 



(35) 

(36) 
(37) 
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Looking for a solution such that $(n*) = O(lnL) implies that 

n* ~aL-lnL + lnhiL a (38) 

where a is a parameter which in principle can be determined. The leading order of h(L) 
is then 

h(L) ~ e~ aL ~ . (39) 

To calculate the width of the peak we note that 

^V) = -(^+^- aL -«^, (40) 

yielding 

r T i 1/2 

A»«~[— ] . (41) 

The condition p(n*)n* An* = 0(1) yields 

p(n*) ~L" 3 / 2 [lnL] 1/2 , (42) 

and the weight w scales as 1/L. 

The scaling behaviour of the peak in p(n) implies that the size of the condensates, 
(38), is extensive. Also, up to logarithmic corrections, the width of the peak (41) scales 
as the square root of its position. Thus the peak is sharp as in the usual condensation 
scenario associated with (1). The number of condensates is 

Q = , 43 

a 

where following our previous discussion we consider only densities p for which q is an 
integer. We will comment on other densities in the next subsection. 

To locate the position of the dip, n min , in p(n), we balance the first two terms in 
the extremum equation (36) which yields 

n ™ (44) 

and 

1 /hiL\ 6 , , 

p{n mm ) ~ - b — • (45) 

n min \ ^ J 

It is interesting to note that an alternative way of obtaining the scaling behaviour 
of the peak of p{n) in the case (4) is by taking the limit k = aL in the results (22,26,29) 
and (31,32) for the algebraic case (3). 
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4-2. Numerical results 

We start by considering a case with integer (p — p c )/a. The results of numerical 
simulations for a system with hopping rates given by (4), size L = 1000, b = 4, a = 1 
and p = 2.5 > p c = 0.5 are displayed in Figure 9. The probability distribution is given 
in Figure 9(a) where we can clearly see a narrow peak around n ~ L = 1000. The 
simulation results are given by the circles and are compared to a numerical calculation 
of the exact canonical expression (7) represented by the full line. Since the critical 
density is 0.5 we expect the excess density 2.5 — 0.5 = 2 to be distributed among two 
condensates. The occupation of each condensate is expected to be approximately equal 
to L — 1000. A typical snapshot is presented in Figure 9(b) where the occupation of 
the lattice is plotted. 




Figure 9. (a) A logarithmic plot of p(n) in the case of algebraic hopping rates obtained 
from numerical simulation, represented by circles (o), compared to that calculated from 
the exact canonical expression (7), represented by the solid line, (b) Site occupations 
in a snapshot of the system obtained from numerical simulation for u(n) given by (4), 
L = 1000, p = 2.5, b = 4 and a = 1. 



We now proceed to examine a case with non-integer (p — p c ) / a. The results of 
numerical simulations for a system with hopping rates given by (4), size L = 1000, 
b = 4, a = 1 and p = 2.3>p c = 0.5 are displayed in Figure 10. The probability 
distribution is given in Figure 10(a) where we can clearly see a rather broad peak 
around n ~ L = 1000. The simulation results are given by the circles and are compared 
to a numerical calculation of the exact canonical expression (7) represented by the full 
line. Since the critical density is 0.5 we expect the excess density 2.3 — 0.5 = 1.8 to 
be distributed among two condensates. Even though the occupation of each condensate 
fluctuates in time between 0.8 and 1 they spend the majority of the time in either and 
not in transition as can be seen in Figure 8. The resulting distribution has two peaks, 
one centred around n\ ps 0.8L = 800 and the second around n\ ~ L = 1000, which 
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yield the appearance of a single broader peak in Figure 10(a). A typical snapshot is 
presented in Figure 10(b) where the occupation of the lattice is plotted. 




Figure 10. (a) A logarithmic plot of p(n) in the case of algebraic hopping rates 
obtained from numerical simulation, represented by circles (o), compared to that 
calculated from the exact canonical expression (7), represented by the solid line . 
(b) Site occupations in a snapshot of the system obtained from numerical simulation 
for u(n) given by (4), L = 1000, p = 2.3, b = 4 and a = 1. 




Figure 11. The canonical probability distribution calculated numerically using (7) 
for L = 3200, in the case of exponential hopping rates with b = 4 a = 1 and p = 1.5, 
1.75 and 2.75 for the full, dotted and dashed lines respectively. 
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In Figure 11, we compare the numerical solution for the canonical distribution p(n) 
for p = 1.5 (integer (p — p c )/a) with that corresponding to p = 1.75 and p = 2.75 (where 
(p — p c )/a is non-integer) for L = 3200, b = 4 and a — 1. For the above parameters 
with a density of p = 1.5 there is a single condensate corresponding to the single peak 
at aL = 3200. On the other hand for p = 2.75 there are two condensates corresponding 
densities (p — p c — a)L = 800 and aL = 3200, resulting in two distinct peaks in the 
distribution function. This is verified in Figure 11. Clearly, the results for non-integer 
(p — p c )/a do not agree with the grand canonical analysis in which there is a single peak 
corresponding to the multiple condensates. 




10° 10 1 10 2 10 3 10° 10 1 10 2 10 3 

n n 



Figure 12. A comparison of the canonical probability (7) and the grand canonical 
probability (10), represented by the dashed and full lines respectively. The 
distributions were calculated numerically for exponential hopping rates with L = 1600, 
b = 4 a = 1 and p = 1.5 and p = 1.75 in figures (a) and (b) respectively. 

Since the steady state properties of the distribution given in Section 4 were 
calculated in the grand canonical ensemble, in which the distribution exhibits only 
a single peak, we wish to compare the distribution obtained in the canonical ensemble 
with the one obtained in the grand canonical ensemble. In Figure 12 the probability 
distribution for the parameters L = 1600, b = 4 and a = 1 was calculated numerically 
both in the grand-canonical ensemble (10) and in the canonical ensemble (7). In 
Figure 12(a) the distributions are compared for p = 1.5, a density in which a single 
condensate exits. In Figure 12(b) the distributions are compared for p = 1.75 a density in 
which two condensates exist. It can be seen in Figure 12(a) and(b) that the distributions 
obtained form the two ensembles are not equivalent; they differ mostly at large n. On 
the other hand, for small n, namely in the fluid phase, the distributions obtained from 
the two ensembles agree. 
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For large occupations not only does the size dependent constraint on the number 
of particles affect the distribution but so does the size dependent term in the hopping 
rate (4). This endogenous size dependence yields a different form for the distribution at 
large n. The height of the peak p(n*) is the same in both ensembles for the two cases 
in Figures 12(a) and (b). However, the dip p(n m i n ) is different for the two ensembles 
even for the single condensate case in Figures 12(a). The scaling behaviour in the 
canonical ensemble of the distribution as a function of the system size can be calculated 
numerically using (7) and compared with the prediction in the grand-canonical ensemble 
(41,42,44,45). 




10 1 10 2 10 3 10 4 

L 

Figure 13. The scaling of the condensate peak p(n*)/ [InL] 1 ^ 2 as a function of the 
system size L in the case of exponential hopping rates (4) with b = 4, p = 2.5 and 
a = 1. The probability p(n*) for each system size L was calculated numerically using 
(7). The data is compared to a linear slope of —3/2. 

In Figure 13 we check that in the case of integer (p — p c )/a the grand canonical 
analysis yields the leading scaling properties of the peak calculated for the canonical 
distribution in finite size systems. In the Figure, p{n*)/ [lnL] 1 ^ 2 is plotted against L. 
A slope which agrees with the expected value of —3/2 is seen, confirming the L~ 3 / 2 
dependence. However to test the logarithmic corrections more detailed numerical data 
are required. 

For the non integer (p — p c )/a case the grand canonical analysis is not valid 
and as a result, the scaling relations developed using the grand canonical analysis do 
not hold. This was confirmed numerically using (7) to calculate the scaling in the 
canonical ensemble and comparing with the prediction in the grand-canonical ensemble 
(41,42,44,45). 

As noted above in the discussion of Figure 12, the canonical and grand canonical 
distributions appear to differ in the dip region. However it may be that scaling is still 
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Figure 14. The scaling of the dip as a function of the system size L in the case of 
exponential hopping rates (4) with b = 4, p = 2.5 and a = 1. The probability p(n m i n ) 
and location of the dip calculated numerically for each system size L using 

(7). In (a) and (b) n m i n and p(n m i n ) are plotted respectively as a a function of Lj ln(L) 
and are compared to a linear fit of the data points. The data points for the various 
system sizes are compare to linear slopes of 1 and —4 in (a) and (b) respectively. 



correctly predicted by the grand canonical analysis. In Figure 14 we check that in 
the case of integer (p — p c )/a the grand canonical analysis yields the leading scaling 
properties of the dip calculated for the canonical distribution in finite size systems. In 
Figure 14(a) the dip location n m i n is plotted against L. A linear fit for the data produced 
a slope of 1.2 ± 0.2. In the figure the data is compared with a slope of 1 corresponding 
to the scaling relation (44). In Figure 14(b) the dip height p(n min ) is plotted against 
L. A linear fit for the data produced a slope of —4.2 ± 0.2. In the figure the data is 
compared with a slope of —4 corresponding to the scaling relation (45). This is not the 
case for the non integer (p — p c )/a since the dip is not well defined as the number of 
humps may be larger than one. 

We conclude that for the case of integer (p — p c )/a the steady state properties 
of the distribution given in Section 4 computed for the grand canonical ensemble are 
able to describe the properties in the canonical ensemble. However, it seems as if the 
analysis breaks down for the case of non integer (p — p c ) /a where the two descriptions 
are qualitatively different in the description of the condensates. 
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5. Dynamics 

5.1. Condensate creation and evaporation timescales 

We are interested in the dynamics of condensate formation for systems described by the 
hopping rates (3,4). Previous studies of condensate formation and dynamics for the case 
(1) have been carried out [6, 7, 13]. We identify two timescales associated with creation 
and evaporation of (meso) condensates in the system. To identify these timescales we 
note that the size of the (meso) condensates fluctuates around n*. These fluctuations 
are occasionally large enough to cause a condensate to evaporate. The typical time for 
a condensate to exist at a given site before it evaporates is termed the evaporation time 
and is denoted by r e . Once a condensate has evaporated the particles are redistributed 
among the other sites and this can result in formation of another condensate. The 
typical time that a given site exists in the fluid phase before a condensate is created at 
that site is termed the creation time and is denoted by r c . 



In Figure 15 we display the time series of the occupation number of a typical site, 
for the algebraic case (3). One clearly sees sharp transitions between the fluid state 
and condensed state. By averaging over a long time series the evaporation and creation 
times may be evaluated. The results are illustrated in Figure 16. Similar plots for the 
exponential case are illustrated in Figure 18 for an integer number of condensates. 

In order to estimate the evaporation time r e and creation time r c we use the 
Arrhenius law. As we have seen, the probability distribution of the occupation of a 
given site is described by an effective potential $(n) as illustrated in Figure 17. In the 
condensed phase <3>(n) exhibits two valleys. The left valley corresponds to fluid states 
and the right valley centred around n* corresponds to condensate states. The potential 
barrier centred around n min corresponds to the dip region in the probability distribution. 
Using the Arrhenius law, a condensate will form on the site once the occupation of the 



n(t) 




Figure 15. The time series for an occupation of a single site for a lattice of size 
L = 100 in the case of algebraic hopping rates (3) with b — A, p = 4, c = 1 and k = 10. 
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Figure 16. The evaporation time r e and creation time r c as a function of the lattice 
size L for an algebraic hopping rate (3) calculated numerically from simulations. The 
result is plotted for algebraic hopping rates with b = 4, p = 3 and k = 5. The results 
were averaged over 100 repetitions for r e and r c in figures (a) and (b) respectively. 
The results are compared with a linear slope (b — l)k/(k + 1) = 2.5 in figure (a) and 
bk/(k + l) = 3.37 in figure (b). 



site crosses the barrier from the left to the right valley. The timescale for creation of a 
condensate is then proportional to the ratio of the probability of being in the fluid to 
the probability of being at the creation threshold n m i n . Thus, since the probability of 
being in the fluid is of the order of 1, we obtain 

T < ~ rih ■ (46) 

For the evaporation time of a condensate, the timescale is proportional to the ratio 
of the probability for a condensate, which is of the order of w, to the probability of 
being at the evaporation threshold which is again at n min , 

(47) 



Note that on general grounds the ratio of the two timescales has to satisfy 

T - ~ J?- (48) 

T c l-W 

where w/(l—w) is the ratio of the probabilities of a site being in a condensate or a fluid 
state. This follows from a steady-state condition that the rate per site of condensate 
creation (1 — w)/r c is balanced by the rate per site of condensate evaporation w/r e . For 
small w, (47) and (46) satisfy (48). 

5.2. Numerical results for algebraically increasing case 

For the algebraic hopping rate (3), the creation and evaporation times scale with the 
system size as 

/ T \ bk/(k+l) 

r c ~ (^) , r e ~ L(»"i)*/(*+i) [lnL] -(i+^)/(*+D . (49) 
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Figure 17. The effective potential $(n) for a lattice of size L = 1600 in the case of 
algebraic hopping rates with 6 = 4, p = 4, c = 1 and k ~ 10. <I>(n) was calculated 
numerically from the exact expression (7) for the hopping rate u{n) given by (3) . 



Using numerical simulations for for b = 4, p = 3 and k = 5 the evaporation time r e 
and creation time r c are plotted in Figure 16 as a function of the lattice size L . In 
Figure 16(a) the evaporation time r e as a function of lattice size is compared on a double 
logarithmic scale to a linear slope of (b — l)k/ (k+ 1) = 2.5. In Figure 16(b) the creation 
time function of lattice size L is compared on a double logarithmic scale to a 

linear slope of bkj (k + 1) = 3.37. 

5.5. Numerical results for exponentially increasing case 

In the following, we separate the discussion of the case of an integer (p — p c )/a from 
that of a non integer (p — p c ) / a. To carry out an analogous analysis for the exponential 
hopping rate (4) with an integer (p — p c )/a. we require an expression for how n m i n 
scales with L. As we have discussed, the hopping rate (4) is the same as the usual 
case (1) but with the introduction of a sharp cut-off at n — aL. Thus we expect p{n) 
to remain the same as that for (1) up to the condensate peak which takes place at 
n = O(L). Therefore n m i n = O(L), as in the usual case (1), and p{n miv ) ~ 1/L b . This 
estimate inserted into (46,47) yields the following leading behaviour for the creation and 
evaporation timescales in the exponential hopping rate case 

r c ~L b , r e ~ L b ~ l . (50) 

We test the leading L behaviour using numerical simulations for b = 3, p = 4 
and a — 1, which corresponds to an integer (p — p c )/a = 3, the evaporation time r e 
and creation time r c are plotted in Figure 18 as a function of the lattice size L . In 
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Figure 18. The evaporation time r e and creation time r c as a function of the lattice 
size L for an exponential hopping rate (4) were calculated numerically from simulations. 
The results are plotted for b = 3, p = 4 and a = 1. The results where averaged over 100 
repetitions for r e and r c in figures (a) and (b) respectively. The results are compared 
to a linear slope b — 1 = 2 in figure (a) and b = 3 in figure (b) . 



Figure 18(a) the evaporation time r e as a function of lattice size is compared on a double 
logarithmic scale to a linear slope of (b — 1) = 2. In Figure 18(b) the creation time r c 
as a function of lattice size L is compared on a double logarithmic to a linear slope of 
6 = 3. These results are consistent with the leading L behaviour derived for the grand 
canonical ensemble (50), although more extensive simulations would be desirable. 

On the other hand, for the case of non-integer (p—p c )/a the the leading L behaviour 
of the creation time r c and of the evaporation time r e which were calculated within 
the grand canonical ensemble are invalid. The Arrhenius law approach is based on a 
distribution with a single peak and a single dip which is not the case for the non-integer 
(p — p c )/a as was shown above. Thus, for the non-integer (p — p c )/a not only is the 
distribution obtained within the grand canonical ensemble not valid, so is the framework 
in which we calculate the evaporation and creation times. 

6. Conclusions 

In this work we have shown how the introduction of non-monotonic hopping rates u(n) 
into the ZRP can produce a condensation transition into many condensates. The 
algebraic choice of hopping rate (3) results in mesocondensates whose number grows 
subextensively with system size L. On the other hand, the exponential choice of hopping 
rate (4) can result in a finite number of extensive condensates. 

Related condensation phenomena, which result in a large number of mesoconden- 
sates, have previously discussed in the context of non-conserving ZRP [9]. Thus our 
results imply that the steady states of such models can be described by a conserving 
ZRP with an effective hopping rate which is non-monotonic. 
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We have analysed the single-site distribution p(n) within the grand canonical 
ensemble and in particular analysed the scaling behaviour of the condensate peak. We 
have computed the scaling behaviour of the position, height, width of the peak and also 
the position and height of the dip between the condensate peak and the part of the 
distribution representing the fluid as illustrated in Figure 2. 

In the case of algebraically increasing hopping rates (3) the predictions of the 
grand canonical analysis are well borne out by numerical computation within the 
canonical ensemble. Ignoring logarithmic corrections the scaling of the condensate peak 
is as follows: the position of the peak scales as n* ~ L k /( k + 1 )-^ the width scales as 
An* ~ L k ^ k+V> ] the height of the peak scales as p(n*) ~ L~ 2k /( k + 1 )^ anc l the weight 
of the peak scales as w ~ [ J ' k /( k + 1 ), The behaviour of w implies that the number of 
meso condensates scales as L 1 /^ 4 " 1 -*. The peak is not sharp since An*/n* — > only as 
a power of InL. Thus we term it a 'weak peak'. The weak peak is similar to that 
determined in the analysis of a non-conserving ZRP [9]. 

For the exponential case (4) the situation is more subtle. For the case of integer 
values of (p — p c )/a, our numerical results suggest that the scaling predicted by the grand 
canonical analysis is correct even though the grand canonical and canonical distributions 
do not appear to coincide. The scaling of the condensate peak, excluding any logarithmic 
corrections, is as follows: the position of the peak scales as n* ~ L; the width scales as 
An* ~ L 1 / 2 ; the height of the peak scales as p(n*) ~ L~ 3 ^ 2 , and the weight of the peak 
scales as w ~ L -1 . The behaviour of w implies that the number of mesocondensates is 
finite. In the exponential case the peak is 'sharp' since An* /n* — > as a power of L. 

For the case of non-integer values of (p — p c )/a, on the other hand, it appears that 
the results of the grand canonical analysis break down. 

It is of interest to compare the scaling form we have determined of the peak in p(n) 
corresponding to (meso) condensates to the usual scenario of condensation into a single 
site exhibited, for example, by the hopping rate (1). For the single condensate case, 
in [10, 11] the condensate peak, denoted p CO nd(n) in that work, has been computed. 
The scaling form of the peak is as follows: the position scales as n* ~ L; the width 
scales as An* ~ L l l 2 for b > 3 and An* ~ L 1 /^- 1 ) for 3 > 6 > 2; the height scales as 
p(n*) ~ L~ 3/2 for b > 3 and p(n*) ~ L~ 6 /(f>-i) for 3 > b > 2; and the weight of the 
peak scales as w ~ L -1 . The peak is sharp, as in the case we have studied of a finite 
number of condensates. However, in the case of a finite number of condensates, our 
grand canonical predictions for the scaling of the width and height of the peak do not 
depend on b. 

We note that for the usual single condensate scenario the results quoted above for 
Pcond(n) were computed within the canonical ensemble [10, 11]. It would be of interest 
to see if our results can be recovered in a canonical calculation. Also, an analysis within 
the canonical ensemble is required to obtain the scaling form the condensate peak in 
the exponential case when (p — p c )/a is non-integer. 

In studying the dynamics of the condensates we identified two timescales, one for 
the creation of a condensate at a given site and the other for the evaporation of a 
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condensate at a given site. The scaling of these timescales with the system size has 
been studied within the phenomenological Arrhenius approach; how these timescales 
scale with system size is determined by the scaling of height of the dip and the weight of 
the condensate peak. We have found that the predictions compare well with numerical 
simulations in the algebraic case. In the exponential case the numerical results are 
consistent with the Arrhenius predictions, but more extensive simulations would be 
desirable. 
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